clc,clear;
close all;
% 种植面积
data=xlsread('2023种植面积data.xlsx');
data(isnan(data))=0;
% 预测销售量
sale=readmatrix('预测销售量.xlsx');
sale=sale(:,3);
sale0=sale;
% 亩产量
chag=readmatrix('亩产量数据表格.xlsx');
chag=chag(:,3:end);
chag0=chag;
% 销售单价
r=readmatrix('销售单价数据.xlsx');
r=r(:,3:end);
r0=r;
% 种植成本
p=readmatrix('种植成本数据.xlsx');
p=p(:,3:end);
p0=p;
% 地区最大面积
S=xlsread('附件1.xlsx'); 
S=S([1:54,27:54],1);
% 最低种植面积
ms=readmatrix('最低种植面积.xlsx');
ms=ms(:,3:end);
%参数
y=7;%年份
m=82;%地区类型
n=41;%作物类型
maxi=100;%最大迭代次数
alap=0.5*0.25+0.5*0.4;%损耗率
T0=300;%初始温度
FT=0.01;%最终温度
a=0.99;%降温率
energy = 1;        %初始能量
energy_d=2;        %能量阈值
ros1=0.5844;
ros2=0.7822;
h=1;
x0=data;%初始数据
bestx=x0;
f0=profit(x0,chag,sale,p,r,alap);
bestf=f0;
alap0=alap;
pop=p0;
ror=r0;
sos=sale0;
for k=1:y
    T=T0; 
    x=x0;
    count=1;
    %进行优化模型求解(二次模拟退火)
    while T>FT
        for i=1:maxi
            %蒙特卡洛随机分布
            lb=1.25;
            ub=1.4;
            alap=(lb+(ub-lb)*rand)*alap;
            for j=1:n
                lb=0.9;
                ub=1.1;
                chag(:,j)=chag(:,j)*(lb+(ub-lb)*rand);
                lb=1.04;
                ub=1.06;
                p(:,j)=p(:,j)*(lb+(ub-lb)*rand);
                if j>=17&&j<=37
                    lb=1;
                    ub=1.05;
                    r(:,j)=r(:,j)*(lb+(ub-lb)*rand);
                elseif j>=38&&j<=40
                    lb=0.95;
                    ub=0.99;
                    r(:,j)=r(:,j)*(lb+(ub-lb)*rand);
                elseif j==41
                    lb=0.94;
                    ub=0.96;
                    r(:,j)=r(:,j)*(lb+(ub-lb)*rand);
                end
            end
            f=profit(x,chag,sale,p,r,alap);
            
            %产出新解
            newx=x;
            %每次对于一类地区中的作物进行一定的扰动
            %(1-26)地区可(1-15)作物,(27-54)可(17-34),产物16只能27-34
            %(55-62)可35-37,(63-78)可38-41,(79-82)可(17-34)
            %分为6类
            op=randi([1,6]);%随机的类
            if op==1
                i0=[1,26];j0=[1,15];
            elseif op==2
                i0=[27,54];j0=[17,34];
            elseif op==3
                i0=[27,34];j0=16;
            elseif op==4
                i0=[55,62];j0=[35,37];
            elseif op==5
                i0=[63,78];j0=[38,41];
            else
                i0=[79,82];j0=[17,34];
            end
       pl=0.5;
      if pl>rand
            if j0==16
                io=randi(i0);
                newx(io,j0)=rand*newx(io,j0)+ms(io,j0);
            else
                io=randi(i0);
                jo=randi(j0,3,1);
                for lk=1:3
                newx(io,jo(lk))=rand*newx(io,jo(lk))+ms(io,jo(lk));
                end
            end
      else
             %交换作物
             if op~=3
                 for hi=min(j0):max(j0)
                     for ki=min(i0):max(i0)
                          if newx(ki,hi)>0
                              msk=newx(ki,hi);
                          else
                              continue;
                          end
                          jo=randi(j0);
                          while jo==hi
                              jo=randi(j0);
                          end
                          mst=newx(ki,jo);
                          newx(ki,jo)=msk;
                          newx(ki,hi)=mst;
                     end
                 end
             end
       end
     %对于豆类提供更容易的概率选中
        po=0.5;
        if po>rand
            jj=1:5;
            ji=[1,5];
            for ki=1:26
                if sum(newx(ki,jj))==0
                    s=randi(ji);
                    newx(ki,s)=ms(ki,s);
                    for kj=6:n
                        if newx(ki,kj)>0
                            newx(ki,kj)=newx(ki,kj)-ms(ki,s);
                            break;
                        end
                    end
                end
            end
            ji=[17,19];
            jj=17:19;
            for ki=27:50
                if sum(newx(ki,jj))==0
                    s=randi(ji);
                    newx(ki,s)=ms(ki,s);
                    for kj=[1:16,20:n]
                        if newx(ki,kj)>0
                            newx(ki,kj)=newx(ki,kj)-ms(ki,s);
                            break;
                        end
                    end
                end
            end
            for ki=51:54
                if sum(newx(ki,jj))==0||sum(newx(ki+28,jj))
                    s=randi(ji);
                    if rand>0.5
                        newx(ki,s)=ms(ki,s);
                    else
                        newx(ki+28,s)=ms(ki+28,s);
                    end
                     for kj=[1:16,20:n]
                        if newx(ki,kj)>0
                            newx(ki,kj)=newx(ki,kj)-ms(ki,s);
                            break;
                        end
                    end
                end
            end
        end
        % 有豆类则增产
        for kops=1:82 
            for jops=[1:5,17:19]
                if newx(kops,jops)>0
                    oplst=[6:16,20:41];
                    chag(kops,oplst)=1.05*chag(kops,oplst);
                    break 
                end
            end
        end
        %关系更新预测销售量
        rs=0;
        ps=0;
        count1=0;
        for ios=1:41
            for jos=1:82
                if r(jos,ios)>0
                    rs=rs+(1-r(jos,ios)/ror(jos,ios));
                    count1=count1+1;
                end
            end
                rs=rs/count1;
                count1=0;
            for jos=1:82
                if p(jos,ios)>0
                    ps=ps+(1-p(jos,ios)/pop(jos,ios));
                    count1=count1+1;
                end
             end
                ps=ps/count1;
                sale(ios)=(1+ros2*ps)*sale(ios);
                sale(ios)=(1-ros1*rs)*sale(ios);
        end
        %约束
        for iops=1:41
            if iops==6||iops==7
                lp=sale(iops)/sos(iops);
                if lp>1.1||lp<1.05
                    sale(iops)=sos(iops)*(rand*0.05+1.05);
                end
            else
                lp=sale(iops)/sale0(iops);
               if lp<0.95||lp>1.05
                    sale(iops)=sale0(iops)*(rand*0.1+0.95);
               end
            end
        end
           %检查是否满足约束条件
        if k==1
        [newx,flag]=yuesu(newx,S,ms,k,x0);
        elseif k==2
           [newx,flag]=yuesu(newx,S,ms,k,x0,y1);
        elseif k==3
            [newx,flag]=yuesu(newx,S,ms,k,y1,y2);
        elseif k==4
            [newx,flag]=yuesu(newx,S,ms,k,y2,y3);
        elseif k==5
            [newx,flag]=yuesu(newx,S,ms,k,y3,y4);
        elseif k==6
            [newx,flag]=yuesu(newx,S,ms,k,y4,y5);
        elseif k==7
            [newx,flag]=yuesu(newx,S,ms,k,y5,y6);  
        end
        if flag==0
            r=r0;
            p=p0;
            chag=chag0;
            sale=sale0;
            alap=alap0;
         continue;
        elseif flag==1&&count==1
            %计算适应值并更新最优值
            newf=profit(newx,chag,sale,p,r,alap);
            %设定初始条件
            bestx=newx;
            bestf=newf;
            f=newf;
            x=newx;
            count=0;
             r1=r;
             p1=p;
             chag1=chag;
             sale1=sale;
             alap1=alap;
        end
            %计算适应值并更新最优值
            newf=profit(newx,chag,sale,p,r,alap);
            delta=newf-f;
            if delta>0||exp(-delta/T)<rand
                x=newx;
                f=newf;
                % 更新全局最优
                if newf>bestf
                     bestx=x;
                     bestf=f;
                     r1=r;
                     p1=p;
                     chag1=chag;
                     sale1=sale;
                     alap1=alap;
                end
            end
            r=r0;
            p=p0;
            chag=chag0;
            sale=sale0;
            % 更新能量
            if delta > 0
                energy = energy + 1;
            end
            alap=alap0;
        end
        % 跳跃式降温
        if energy > energy_d
            T = T * (a ^ (energy_d*sqrt(energy-energy_d)));  % 快速降温
        else
            T = T * a;  % 正常降温
        end
        % 重置能量
           energy = 1;
           h=h+1;
    end
% 二重退火：重启退火过程
    T=T0;
    energy = 0;
    r=r0;
    p=p0;
    chag=chag0;
    sale=sale0;
    count=1;
    while T>FT
        for i=1:maxi
            %蒙特卡洛随机分布
            lb=1.09;
            ub=1.11;
            alap=(lb+(ub-lb)*rand)*alap;
            for j=1:n 
                lb=0.9;
                ub=1.1;
                chag(:,j)=chag(:,j)*(lb+(ub-lb)*rand);
                lb=1.04;
                ub=1.06;
                p(:,j)=p(:,j)*(lb+(ub-lb)*rand);
                if j>=17&&j<=37
                    lb=1;
                    ub=1.05;
                    r(:,j)=r(:,j)*(lb+(ub-lb)*rand);
                elseif j>=38&&j<=40
                    lb=0.95;
                    ub=0.99;
                    r(:,j)=r(:,j)*(lb+(ub-lb)*rand);
                elseif j==41
                    lb=0.94;
                    ub=0.96;
                    r(:,j)=r(:,j)*(lb+(ub-lb)*rand);
                end
            end
            f=profit(x,chag,sale,p,r,alap);
            
            %产出新解
            newx=x;
            %每次对于一类地区中的作物进行一定的扰动
            %(1-26)地区可(1-15)作物,(27-54)可(17-34),产物16只能27-34
            %(55-62)可35-37,(63-78)可38-41,(79-82)可(17-34)
            %分为6类
            op=randi([1,6]);%随机的类
            if op==1
                i0=[1,26];j0=[1,15];
            elseif op==2
                i0=[27,54];j0=[17,34];
            elseif op==3
                i0=[27,34];j0=16;
            elseif op==4
                i0=[55,62];j0=[35,37];
            elseif op==5
                i0=[63,78];j0=[38,41];
            else
                i0=[79,82];j0=[17,34];
            end
       pl=0.5;
      if pl>rand
            if j0==16
                io=randi(i0);
                newx(io,j0)=rand*newx(io,j0)+ms(io,j0);
            else
                io=randi(i0);
                jo=randi(j0,3,1);
                for lk=1:3
                newx(io,jo(lk))=rand*newx(io,jo(lk))+ms(io,jo(lk));
                end
            end
      else
             %交换作物
             if op~=3
                 for hi=min(j0):max(j0)
                     for ki=min(i0):max(i0)
                          if newx(ki,hi)>0
                              msk=newx(ki,hi);
                          else
                              continue;
                          end
                          jo=randi(j0);
                          while jo==hi
                              jo=randi(j0);
                          end
                          mst=newx(ki,jo);
                          newx(ki,jo)=msk;
                          newx(ki,hi)=mst;
                     end
                 end
             end
       end
     %对于豆类提供更容易的概率选中
        po=0.5;
        if po>rand
            jj=1:5;
            ji=[1,5];
            for ki=1:26
                if sum(newx(ki,jj))==0
                    s=randi(ji);
                    newx(ki,s)=ms(ki,s);
                    for kj=6:n
                        if newx(ki,kj)>0
                            newx(ki,kj)=newx(ki,kj)-ms(ki,s);
                            break;
                        end
                    end
                end
            end
            ji=[17,19];
            jj=17:19;
            for ki=27:50
                if sum(newx(ki,jj))==0
                    s=randi(ji);
                    newx(ki,s)=ms(ki,s);
                    for kj=[1:16,20:n]
                        if newx(ki,kj)>0
                            newx(ki,kj)=newx(ki,kj)-ms(ki,s);
                            break;
                        end
                    end
                end
            end
            for ki=51:54
                if sum(newx(ki,jj))==0||sum(newx(ki+28,jj))
                    s=randi(ji);
                    if rand>0.5
                        newx(ki,s)=ms(ki,s);
                    else
                        newx(ki+28,s)=ms(ki+28,s);
                    end
                     for kj=[1:16,20:n]
                        if newx(ki,kj)>0
                            newx(ki,kj)=newx(ki,kj)-ms(ki,s);
                            break;
                        end
                    end
                end
            end
        end
       % 有豆类则增产
        for kops=1:82 
            for jops=[1:5,17:19]
                if newx(kops,jops)>0
                    oplst=[6:16,20:41];
                    chag(kops,oplst)=1.05*chag(kops,oplst);
                    break 
                end
            end
        end
        %关系更新预测销售量
        rs=0;
        ps=0;
        count1=0;
        for ios=1:41
            for jos=1:82
                if r(jos,ios)>0
                    rs=rs+(1-r(jos,ios)/ror(jos,ios));
                    count1=count1+1;
                end
            end
                rs=rs/count1;
                count1=0;
            for jos=1:82
                if p(jos,ios)>0
                    ps=ps+(1-p(jos,ios)/pop(jos,ios));
                    count1=count1+1;
                end
             end
                ps=ps/count1;
                sale(ios)=(1+ros2*ps)*sale(ios);
                sale(ios)=(1-ros1*rs)*sale(ios);
        end
        %约束
        for iops=1:41
            if iops==6||iops==7
                lp=sale(iops)/sos(iops);
                if lp>1.1||lp<1.05
                    sale(iops)=sos(iops)*(rand*0.05+1.05);
                end
            else
                lp=sale(iops)/sale0(iops);
               if lp<0.95||lp>1.05
                    sale(iops)=sale0(iops)*(rand*0.1+0.95);
               end
            end
        end
        %检查是否满足约束条件
        if k==1
        [newx,flag]=yuesu(newx,S,ms,k,x0);
        elseif k==2
           [newx,flag]=yuesu(newx,S,ms,k,x0,y1);
        elseif k==3
            [newx,flag]=yuesu(newx,S,ms,k,y1,y2);
        elseif k==4
            [newx,flag]=yuesu(newx,S,ms,k,y2,y3);
        elseif k==5
            [newx,flag]=yuesu(newx,S,ms,k,y3,y4);
        elseif k==6
            [newx,flag]=yuesu(newx,S,ms,k,y4,y5);
        elseif k==7
            [newx,flag]=yuesu(newx,S,ms,k,y5,y6);  
        end
        if flag==0
            r=r0;
            p=p0;
            chag=chag0;
            sale=sale0;
            alap=alap0;
         continue;
        elseif flag==1&&count==1
            %计算适应值并更新最优值
            newf=profit(newx,chag,sale,p,r,alap);
            %设定初始条件
            bestx=newx;
            bestf=newf;
            f=newf;
            x=newx;
            count=0;
            r1=r;
            p1=p;
            chag1=chag;
            sale1=sale;
            alap1=alap;
        end
            %计算适应值并更新最优值
            newf=profit(newx,chag,sale,p,r,alap);
            delta=newf-f;
            if delta>0||exp(-delta/T)<rand
                x=newx;
                f=newf;
                 % 更新全局最优
                if newf>bestf
                     bestx=x;
                     bestf=f;
                     r1=r;
                     p1=p;
                     chag1=chag;
                     sale1=sale;
                     alap1=alap;
                end
            end
            r=r0;
            p=p0;
            chag=chag0;
            sale=sale0;
            % 更新能量
            if delta > 0
                energy = energy + 1;
            end
            alap=alap0;
        end
        % 跳跃式降温
        if energy > energy_d
            T = T * (a ^ (energy_d*sqrt(energy-energy_d)));  % 快速降温
        else
            T = T * a;  % 正常降温
        end
        % 重置能量
           energy = 1;
           h=h+1;
    end
    disp(bestf)
    if k==1
        y1=bestx;
    elseif k==2
        y2=bestx;
    elseif k==3
        y3=bestx;
    elseif k==4
        y4=bestx;
    elseif k==5
        y5=bestx;
    elseif k==6
        y6=bestx;
    elseif k==7
        y7=bestx;
    end
    bestx=x0;
    bestf=f;
    f=f0;
    r=r1;
    p=p1;
    chag=chag1;
    sale=sale1;
    alap=alap1;
    pop=p1;
    ror=r1;
    sos=sale1;
end
function f=profit(x,ch,y,p,r,alpha)
   f=0;
   for i=1:41
        if sum(x(:,i).*ch(:,i))-y(i)>0
            t=1;
        else
            t=0.5;
        end
        f=f+sum((1-alpha)*x(:,i).*ch(:,i).*r(:,i)-x(:,i).*p(:,i))-sum(t*(1-alpha)*r(:,i).*(sum(x(:,i).*ch(:,i))-y(i)));
   end
end

function [x,flag]=yuesu(x,S,ms,k,yx1,yx2)
   [m,n]=size(x);
   flag=1;
   for i=1:m
           %种植面积约束
           if x(i,:)>S(i)
               flag=0;
               return
           end
       for j=1:n
           %正向约束
           if x(i,j)<0
               flag=0;
               return
           end
       end
   end 
   %最小种植面积与地区约束
   for i=1:26  %第一大类
      for j=1:n
          if j>15
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif j==11&&x(i,j)>0
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
           if  j<=15
              if k==1
                  if yx1(i,j)*x(i,j)>0
                    x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                      x(i,j)=0.9*x(i,j);
                  end  
              end
           end
          if j<=15&&j~=11&&x(i,j)>0
              yo=[1:10,12:15];
              opl=1;
              to1=zeros(1,length(yo));
             for ki=yo
                 if x(i,ki)>0
                     to1(opl)=ki;
                     opl=opl+1;
                 end
             end
             if sum(x(i,to1(1:opl-1)))<sum(ms(i,to1(1:opl-1)))
                    flag=0;
                    return
              end
          end 
      end
   end
   for i=27:34  %第二大类
      for j=1:n
          if j<16||j>34
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif x(i,j)>0&& (j==16||j==23||j==24||j==25||j==26||j==30||j==31||j==32||j==33||j==34)
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
          if  j>=16&&j<=34
              if k==1
                  if yx1(i,j)*x(i,j)>0
                    x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                      x(i,j)=0.9*x(i,j);
                  end  
              end
          end
          if (j>=27&&j<=29)||(j>=17&&j<=22)&&x(i,j)>0
              yo=[17:22,27:29];
              opl=1;
             to2=zeros(1,length(yo));
             for ki=yo
                 if x(i,ki)>0
                     to2(opl)=ki;
                     opl=opl+1;
                 end
             end
             if sum(x(i,to2(1:opl-1)))<sum(ms(i,to2(1:opl-1)))
                    flag=0;
                    return
              end
         end
      end
   end
   for i=35:54  %第三大类
      for j=1:n
          if j<17&&j>34
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif x(i,j)>0&&(j==23||j==24||j==25||j==26||j==30||j==31||j==32||j==33||j==34)
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
           if  j>=17&&j<=34
              if k==1
                  if yx1(i,j)*x(i,j)>0
                       x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                       x(i,j)=0.9*x(i,j);
                  end  
              end
           end
            if ((j>=27&&j<=29)||(j>=17&&j<=22))&&x(i,j)>0
               yo=[17:22,27:29];
              opl=1;
              to3=zeros(1,length(yo));
             for ki=yo
                 if x(i,ki)>0
                     to3(opl)=ki;
                     opl=opl+1;
                 end
             end
             if sum(x(i,to3(1:opl-1)))<sum(ms(i,to3(1:opl-1)))
                    flag=0;
                    return
              end
           end
      end
   end
   for i=55:62  %第四大类
      for j=1:n
          if j<35&&j>37
              if x(i,j)>0
                  flag=0;
                  return
              end
          end
          if  j>=35&&j<=37
              if k==1
                  if yx1(i,j)*x(i,j)>0
                       x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                       x(i,j)=0.9*x(i,j);
                  end  
              end
          end
          if  j>=35&&j<=37&&x(i,j)>0
              yo=35:37;
              opl=1;
             to4=zeros(1,length(yo));
             for ki=yo
                 if x(i,ki)>0
                     to4(opl)=ki;
                     opl=opl+1;
                 end
             end
             if sum(x(i,to4(1:opl-1)))<sum(ms(i,to4(1:opl-1)))
                    flag=0;
                    return
             end
          end
      end
   end
   for i=63:78  %第五大类
      for j=1:n
          if j<38&&j>41
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif x(i,j)>0&&  j==41
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
           if j>=38&&j<=41
              if k==1
                  if yx1(i,j)*x(i,j)>0
                      x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                         x(i,j)=0.9*x(i,j);
                  end  
              end
           end
            if j>=38&&j<=40&&x(i,j)>0
              yo=38:40;
             opl=1;
             to5=zeros(1,length(yo));
             for ki=yo
                 if x(i,ki)>0
                     to5(opl)=ki;
                     opl=opl+1;
                 end
             end
             if sum(x(i,to5(1:opl-1)))<sum(ms(i,to5(1:opl-1)))
                    flag=0;
                    return
              end
           end
      end
   end
   for i=79:82  %第六大类
      for j=1:n
          if j<17&&j>34
              if x(i,j)>0
                  flag=0;
                  return
              end
          elseif x(i,j)>0&& (j==23||j==24||j==25||j==26||j==30||j==31||j==32||j==33||j==34)
              if x(i,j)<ms(i,j)
                 flag=0;
                 return
              end
          end
          if  j>=17&&j<=34
              if k==1
                  if yx1(i,j)*x(i,j)>0
                      x(i,j)=0.9*x(i,j);
                  end
              else
                  if yx2(i,j)*x(i,j)>0
                        x(i,j)=0.9*x(i,j);
                  end  
              end
          end
          if ((j>=27&&j<=29)||(j>=17&&j<=22))&&x(i,j)>0
               yo=[17:22,27:29];
              opl=1;
             to6=zeros(1,length(yo));
             for ki=yo
                 if x(i,ki)>0
                     to6(opl)=ki;
                     opl=opl+1;
                 end
             end
             if sum(x(i,to6(1:opl-1)))<sum(ms(i,to6(1:opl-1)))
                    flag=0;
                    return
              end
          end
      end
   end
   %三年豆类约束
   j=1:5;
   for i=1:26
       if k>1
           if sum(yx1(i,j))+sum(yx2(i,j))+sum(x(i,j))==0
               flag=0;
               return
           end
       end
   end
   j=17:19;
   for i=27:50
       if k>1
           if sum(yx1(i,j))+sum(yx2(i,j))+sum(x(i,j))==0
               flag=0;
               return
           end
       end
   end
   for i=51:54
       if k>1
           if sum(yx1(i,j))+sum(yx2(i,j))+sum(x(i,j))==0&&sum(yx1(i+28,j))+sum(yx2(i+28,j))+sum(x(i+28,j))==0
               flag=0;
               return
           end
       end
   end
end


